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ABSTRACT 

The  range  dependent  boundary  condition  imposed  by  a rough  ice 
surface  serves  to  couple  the  normal  modes  of  an  underwater  a- 
coustic  field.  Beginning  with  the  parabolic  wave  equation*  the 
solution  is  expanded  in  terms  of  normal-mode  depth  functions 
modulated  by  random  amplitudes.  A system  of  master  equations 
for  the  quantities  <an  (r)  an*  (r)  > is  derived*  where  the  a^r) 
are  the  normal-mode  amplitudes,  and  the  brackets  indicate  an  en- 
semble average.  The  master  equations  determine  the  mean  power 
in  each  mode  as  a function  of  range*  describing  the  transfer  of 
energy  between  modes.  Explicit  relations  for  the  coupling  coef- 
ficients are  obtained  in  terms  of  the  spectrum  of  the  rough  sur- 
face. Transmission  loss*  intensity  redistribution  and  the  ap- 
proach to  equipartition  of  energy  as  functions  of  range  and  depth 
are  obtained  numerically. 


■BDHMBBi 


■ 
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I.  INTRODUCTION 


The  problem  of  acoustic  scattering  from  the  rough  ice  cover 


of  the  Arctic  Ocean  has  been  of  considerable  interest  for  a long 


time  because  the  effect  is  unavoidable.  To  the  present,  the  mul- 


tiple scattering  phenomenon  has  not  been  considered  with  the  ex- 


ception of  recent  work  by  Kryazhev,et  alA.  The  present  work  will 


be  compared  in  detail  with  the  Soviet  work  in  a future  report. 


The  present  calculation  makes  use  of  the  method  of  parabolic 


equations  (forward  scattering,  long  range,  kQr>>l)  introduced  to 


underwater  acoustics  by  P.  D.  Tappert*  to  study  long-range  propa- 


gation. The  medium  we  consider  is  the  Arctic  Ocean  whose  sound 


speed  profile  is  taken  to  be  bilinear  . The  shape  of  the  inter- 


face between  the  under-ice  cover  and  the  water  is  n(r),  where  n (r) 


is  a random  ice  roughness  function  represented  by  its  wavenumber 


distribution.  The  integral  of  this  wavenumber  distribution  is  the 


r.m.s.  under-ice  amplitude  squared.  The  form  of  the  wavenumber 


distribution  is  motivated  by  previous  research  carried  out  by 


Hibler,  et  al  and  is  taken  here  to  be  Gaussian. 


The  rough  surface  causes  the  initial  acoustic  intensity  to  be 


redistributed  as  a function  of  depth.  It  is  also  found  that  the 


energy  per  mode  approaches  equipartition  as  the  observation  range 


approaches  long  ranges. 


The  reader  is  referred  to  Ref.  5 to  obtain  further  details  to 


the  calculations  shown  below.  It  is  noted  here  that  this  is  a pre- 


liminary report  to  be  followed  by  a report  including  comparison 


to  experimental  data. 


Page  3 


II.  Theory. 


Consider  a point  source  at  point  (0,zQ)  in  a cylindrical 
coordinate  system,  emitting  an  acoustic  signal  at  angular  fre- 
quency u.  The  problem  is  to  describe  the  average  effect  of  a 
randomly  distributed  rough  surface  on  the  distribution  of 
energy  within  the  ocean  waveguide.  Of  particular  interest  is 
the  average  intensity  as  a function  of  range  and  depth 
I(r,z)  - <p(r,z)p* (r,z) > where  p(r,z)  is  the  pressure  ampli- 
tude. 

We  begin  with  the  equation  satisfied  by  the  envelope 


modulating  the  cylindrical  plane  wave,  i.e.  the  parabolic  wave 


equation: 


2ik03*(r,z)/3r  + 32i|>(r,z)/3z2 


(1) 


+ k2(n2(z)  - l)ij>  (r,z)  - 0 , 


where  the  envelope  function  tp(r,s)  satisfies  the  following 
boundary  conditions: 


\|i(r,n(r) ) " 0 , 


3t|»(r,z)/3z 


Z"H 


- o , 


(la) 

(lb) 


*(0,z)  - (8ir)c0)"1/2«(2  - 2q)  , 


(lc) 


where  a (r)  represents  the  randomly  distributed  boundary  between 
the  water  and  under-ice  surface,  and  H is  -the  ocean  depth. 


(2) 


Performing  the  coordinate  transformation5 
r'  « r , 

- s - n(r)  , 

and  rewriting  <p(  r,z)  as  iMr'  ,z'  )exp(i8  (r'  ,z' ) ) , Eq.  (1)  be- 
comes (dropping  the  primes,  suppressing  the  arguments  for  con- 
venience and  imposing  conditions  on  6(r',z')) 

i3*/3r  + (l/2k0) 32$/3z2  + (k0/2)(n2(z)  -1)$ 


- kQz32Ti/3r25;  « 0 , (3) 

$ (r,0)  - 0 , (3a) 
3*/3z|z-H  “ 0 ' (3b) 
? <0,z)  - (8irk 0)"1/2d(z  - iQ)  , (3c) 


where  zQ  ■ zQ  - n(0).  The  random  boundary  condition  (Eq.  (la)) 
becomes  a "flat  surface"  boundary  condition  at  the -expense  of 
introducing  an  additional  perturbative  term  in  the  parabolic 
wave  equation. 

The  envelope  function,  $,  is  expanded  in  terms  of  a 
complete  set  of  "states"  or  normal  modes,  defined  by  the  unper- 
turbed sound-speed  profile.  Each  normal  mode  is  modulated  by 
a random  amplitude  indicating  the  effect  of  the  surface  on 
acoustic  propagation: 
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ip  (r,z)  ■ Xa  an(r>  ♦n<*)exp[i(Xn  - kQ/2)r]  , 


where 


t(l/2k0)32/3z2  + (k0/2)n2(z)]«Pn(z)  - *n4>n(z)  ' 


n 

f ♦„<*>♦„ 


(z)dz  - 6__  . 

run 


4 " 


Inserting  Eq.  (4)  into  Eq.  (3)  and  using  the  orthonormality  of 
the  normal -mode  depth  function,  one  finds,  upon  equating  the 
coefficients  of  4>n(z)exp[idn  - kQ/2)r],  that  the  random 
amplitudes  satisfy  the  coupled  mode  equations, 

da  (rVSr  - i \ c (r)exp(il r)a_(r)  , (5) 

n ~ itin  inn  m 

n 


where 


, 

mn  m n 


H 

(r)  - | 4m(z) H(r,z)<|>n 


(z)dz 


with 


- -k0zmn»2,ier>/9r2  ' 


H (r,z)  - -kQz32n(r)/3r2  , 


Z - Z 
mn  nm 


n 


(z)dz  , 


is  the  "dipole  moment".  The  coupling  coefficients,  c^tr) 
are  real  and  symmetric.  The  coupling  coefficients  are 
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expressed  in  terns  of  a Fourier  integral  representation.  Im- 
posing the  conditions  that  the  mean  amplitude  of  the  surface 
roughness  is  zero  and  that  correlations  of  coupling  coeffi- 
cients be  spatially  stationary,  we  obtain  specific  conditions 
on  the  Fourier  amplitudes.  Integrating  Eq.  (5)  from  r to 
r + Lx,  where  Lx  « x,  to  second  order  in  the  coupling  coeffi- 
cients and  anticipating  use  of  the  above  mentioned  statistical 
properties  of  the  Fourier  amplitudes,  an  approximate  expression 

for  the  a (r)  is  obtained.  One  then  constructs  the  average 
n 

correlation  amplitude  between  modes  n and  m applying  the  random 
phase  approximation  (this  quantity  appears  in  the  average 
intensity) : 


A (r)  - <a  (r)a*(r)>  - < | a (r)  1 2>«5__ 
nm  n m n run 


A (r)  6 
n nm 


(6) 


Using  the  expression  for  an(?)  obtained  by  integrating  the 
coupled  mode  equations,  we  obtain  the  master  equations  govern- 
ing the  evolution  of  AR(r)  as  one  steps  out  in  range: 


where 


dAn (r) /dr 


r"  - r8  | z 

nm  mn  0 1 nm 


CjAfr)  - A (r)  ] , 
nm  m n 


(7) 


2 .4 


k P.  (k) 


k«A 


nm 


with  P^(k)  being  the  under-ice  surface  wavenumber  distribution 
whose  integral  over  all  wavenumbers  is  the  r.m.s.  surface 
height  squared. 


~ 8 t 7*7  : ' 


— ; 
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Eq.  (7)  is  i tenable  to  physical  interpretation.  If  mode 
n has  the  largest  energy  content  at  some  range,  then  Eq.  (7) 
tells  us  that  as  we  step  out  in  range,  mode  n will  lose  energy 
to  the  other  modes  of  the  system  (dAR/dr  < 0).  If  in  the 
process  of  releasing  this  energy,  mode  n'  ^ n obtains  a 
larger  share  of  the  energy  than  the  other  modes,  then  a mode 
n'  will  release  energy  to  the  other  modes  (dAn,/dr  < 0).  This 
mixing  of  energy  continues  as  one  marches  out  in  range  until 
sufficient  multiple  scattering  has  occurred  such  that  the 
energy  becomes  equipartiioned.  In  that  case  Am(r)  An(r)  or 
dA  (r)/dr  0.  From  Eq.  (7),  summing  over  all  modes,  one 

® r-HB 

finds 

£ cS  v*n  - 0 ■ 

or 

l A (r)  - constant  = , A (0)  . (8) 

n n 

In  view  of  Eq.  (8)  and  the  interpretation  of  Eq.  (7)  we  con- 
clude 

A (r)  -*■  s-  for  all  n - 1,N  , (8a) 

n r-*» 

where  N is  the  total  number  of  modes  in  the  system  (restricted 
to  include  all  propagating  modes) . 

Eq.  (8)  is  a statement  of  energy  conservation,  which 
should  be  expected  since  we  are  considering  initially  a medium 
with  no  loss  mechanisms. 


III.  CALCULATIONS 


We  begin  with  the  master  equations 


d*  (rj/ta  =■  j OVr)  - Vr>) 

“ m^n 


(9) 


Defining 


run 


t1  - uc  in  r;» 


(10) 


where  £ is  a dummy  index  r Eq.  (9)  can  be  rewritten  as 


dVr,/<Jr  - Z Vr) 

m 


(11) 


or  in  matrix  notation 

✓ 

dA(r)/dr  =*  A A(r)  . 


(11a) 


Defining  P to  be  the  matrix  whose  columns  are  comprised  of  the 
eigenvectors  of  A with  eigenvalues  d,  we  define  the  quantity 


B(r) 


PT  A(r)  . 


(12) 


Employing  the  orthogonality  of  P (PTP  • P PT  » 1) , 
Eq.  (11a)  becomes 

dB(r)/dr  - D B(r) 


(13) 


where  D ■ PT  A 


P is  diagonal,  the  diagonal  elements  being 
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the  eigenvalues  of  A.  If  we  rewrite  Eg.  (13)  as 


dB  (r)/dr  » l d B (r)  (13a) 

m mn  n 

n 

using  - “ldni^mn  (Gerschgorin ' s Theorem6),  each 

element  of  Eq.  (13a)  can  be  solved 


Bm(r)  » exp  Hdmlr)  Bm(0) 


(13b) 


Inverting  the  transformation  Eq.  (12) , Eq.  (13b)  yields 

An(r)  “ E pnk  e3cp  PkA  Ai(0)  * 

k , A 

Calculating  An(r  + R)  we  find 

Vr+R)  “ J " Pnk  exp  (“ldk>R)  PkA  AA(r). 

jC  § 


(14) 


(15) 


Note  that  Eq.  (15)  is  an  exact  relation,  R is  not  limited  in 
magnitude.  Eq.  (15)  describes  how  to  advance  the  average  cor- 
relation amplitude  to  range  r + R knowing  its  value  at  range  r. 

In  particular,  if  we  know  the  initial  value  of  Afl,  i.e.  Aq(0), 
we  can  calculate  A (r)  for. any  range  r > 0 using  Eq.  (15) . 

XX 

Numerically,  in  order  to  avoid  recalculating  a double  summa- 
tion, which  is  time  consuming  and  uses  unnecessary  storage,  we 
define  the  matrix  CnA(R)  by 


c»t(s>  • l pnk  «P  . 


(16) 
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This  matrix  can  be  calculated  once  for  each  roughness  setting  and 
stored.  Eq.  (IS)  is  recast  into  the  form 


Vr+B)  ' ? CnltR>  Vr>  . 


In  Eq.  (17) , C can  be  regarded  as  a "translation  operator1 
The  quantity  of  interest  here  is 


I(r,z)  - < |p(r,z,t) |2>  - J (r,z)/r 


where 


J(r,z)  - l An(r)  <fr*(z)/  Z An(r) 


(18a) 


I(r,z)  is  called  the  average  intensity  and  does  not  depend  on  t 
within  the  context  of  our  model.  Removing  the  cylindrical  spread- 
ing factor,  1/r,  J(r,z)  shows  how  the  rough  surface  affects  the 
propagation  of  the  acoustic  signal.  We  make  use  of  the  initial 


condition  to  find 


a^O)  - (<2irA0>1/2/4ir)  W 


(19a) 


and  thus 


An(0)  - (l/8irk0)  *n2(z0)  . (19b) 


Knowing  this  initial  condition  completely  determines  the  average 
intensity. 

• We  now  discuss  calculation  of  the  normal-mode  depth  functions, 
•n(z).  The  depth  functions  are  calculated  by  a multiple  shooting 
method.  The  computer  code  that  we  use  was  developed  by  L.  B. Dozier ' 
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applied  to  the  ocean  waveguide.  Similar  techniques  were  used 
before  in  solving  one  dimensional  radial  Schroedinger  equations. 
Boundary  conditions  are  imposed  at  z ■ 0,  H,  then  the  normal 
mode  depth  equations  are  integrated  (by  a Numerov  scheme)  from 
the  surface  into  the  deep  ocean  and  simultaneously  from  the 
ocean  floor  upward.  The  integrations  continue  to  a matching 
point  somewhere  in  between.  In  the  process  of  integrating  a trial 
eigenvalue  is  used  (it  is  initially  estimated  using  a WKB  approach) . 
A comparison  of  the  solution  at  the  matching  point  is  performed  to 
properly  match  the  upward  and  downward  integrations.  This  also 
allows  one  to  improve  the  initial  guess  made  for  the  eigenvalue 
for  that  particular  mode.  The  improvement  of  the  eigenvalue  esti- 
mate is  carried  out  by  an  iterative  process  and  continues  until  con- 
vergence to  some  preset  error  is  reached. 

The  diagonalization  of  the  matrix  A and  calculation  of  its 
eigenvalues  are  performed  by  subroutines  developed  by  the  Argonne 
National  Laboratory  in  a program  package  called  EISPACK.  The  sub- 
routines are  called  TRED2  and  IMTQL2.  TRED2  reduces  a real  sym- 
metric matrix  into  a symmetric  tridiagonal  matrix.  IMTQL2  takes 
the  tridiagonal  matrix  and  calculates  its  eigenvalues  and  eigen- 
vectors. It  forms  the  orthogonal  matrix  P whose  columns  are  the 
eigenvectors  of  A. 

The  above  subroutines  (TRED2  and  IMTQL2)  are  employed  in  a 
computer  code  to  perform  the  manipulations  leading  to  the  con- 
struction of  A and  the  subsequent  implementations  of  Eqs.  (17) 
and  (18) . The  results  of  the  calculations  of  the  normal  modes 
are  read  in  from  an  external  tape. 
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IV.  APPLICATION  TO  ARCTIC  ICE  SCATTERING 
A.  MODELS 

We  represent  the  velocity  profile  as  bilinear: 


c(z)/cn  - a + bz  , 0 < z < z 

u — — in 


(20) 


« c + dz  , z < z < H , 
in  — — 


where  a - 1,  b - 3.48765  x lo“5m-1,  c - 1.006106525 

and  d - 1.67741  10~®  m”1,  z Is  the  depth  at  which  the  two  po si- 

in 

tlve  gradients  b and  d meet;  z_  ■ 337.332  m.  H Is  the  Arctic 

In 

Ocean  depth,  taken  to  be  4 kilometers. 

The  multiple  shooting  method  developed  to  calculate  normal 
modes  in  the  non-ice  covered  ocean  has  been,  adapted  to  calculate 
the  normal  modes  for  the  bilinear  profile.  Because  of  the  struc- 
ture of  the  Arctic  velocity  profile,  there  are  only  RSR  modes 
(RSRBR  modes  again  will  not  be  considered) . Therefore  the  scat- 
tering will  cause  more  rapid  energy  transfer  among  the  RSR  modes 
than  in  the  nonpolar  oceans.  Thus  equipartition  is  reached  sooner 
as  compared  with  nonpolar  ocean  scattering. 

An  important  point  to  be  made  is  that  due  to  the  nature  of 
the  bilinear  profile,  it  is  necessary  to  use  a different  initial 
guess  for  the  eigenvalues  from  the  one  used  for  the  canonical  pro- 
file of  the  deep  sound  channel.  To  estimate  the  initial  eigenvalue 
we  use  the  fact  that  the  gradients  are  very  small  for  the  bilinear 
profile  and  approximate  the  profile  as  linear.  We  then  go 
through  the  WKB  analysis  for  this  initial  guess,  using  the  con- 


t 
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nection  formulas,  and  obtain  for  our  guess. 


(0) 


(3bv/2kft)2/3  (2n  + 3/2) 2/3 


(21) 


A -fundamental  difference  between  Arctic  ice  scattering  and  ocean 
surface  scattering  is  the  surface  wavenumber  spectrum.  We  take 
the  under- ice  surface  wavenumber  distribution  (keel  distribution) 
to  be  Gaussian4  with  r.m.s.  wave  height  a and  correlation  length  L: 


Px(k)  - (c2L/Tr1/2)  exp  (-k2L2) 


(22) 


The  value  of  L is  based  on  experimental  observation  made  by 

g 

0.  I.  Diachok.  He  measured  the  order  of  10  ice  ridges  per 
kilometer.  Thus  we  take  L to  be  of  order  100  meters. 

The  corresponding  master  equations  are 


dAn(r)/dr  - I rnn(Am(r)  - An(r)} 

wrn 


(23) 


where  now 


r*  - r 

run 


nm 


-2.(knVu^A1''2]  «cp  (-X2„L2)  |Z,  i2 


nm 


nm1 


(24) 


where  ahd  Z have  the  same  meaning  as  before. 

nut  nm 

Before  going  to  the  results,  we  mention  the  similarities 
and  differences  between  our  formalism  and  that  given  by  Kryaxhev, 
at  al.1  They  assume  a fixed  boundary,  allow  for  volume  attenuation 
in  the  ice  and  bottom,  formally  consider  2-0  transverse  phenomena 


Page  14 

and  their  integral  "transport  equation",  satisfied  by  A (r) , must 
be  calculated  at  every  range  point  r.  In  our  formalism,  the  sur- 
face can  be  moving  or  fixed.  The  moving  surface  application  al- 
lows qs  to  calculate  long-range  power  spectra . We  have  not  yet 
described  dissipative  mechanisms  in  our  formalism,  but  this  has 
in  fact  been  done  as  explained  below.  Our  formalism  is  much  more 
efficient  when  it  cones  to  calculating  An(r).  Finally,  they  cal- 
culate a quantity  called  the  propagation  anomaly  (PA)  related  to 
our  transmission  loss  (TL)  by 

PA  - TL  ♦ 20  log  (r/r^  (25) 

where  r^  ■ 1 yd. 

Attenuation  due  to  surface  scattering  has  been  studied  by  Marsh, 
Schulkin  and  Kneale10.  The  faetor 

sine/r^ 

in  the  attenuation  constant  is  evaluated  for  a limiting  ray  whose 
vertex  depth  is  1 km. 

Since 

sine-ve  - {1-(C  /C(»)  )2}1/2 

O 

and 

C'/Cix)  • (1  ♦ d*)-1 


with  r^  being  determined  by  ray  theory 


Page  15 


(r^  ■ 29Q/d,  eQ  initial  ray  direction) 
it  is  found  that 

sine/^  ■ • d/2  , 

which  is  independent  of  r^.  Thus  an  empirical  relation  ex- 
pressing ice  scattering  per  unit  distance  (dB/km)  is  obtained, 
which  depends  upon  the  source  frequency,  f,  and  the  r.m.s.  sur- 
face roughness,  a s 

SMSK  * 9-3373  x 10”6  f3/2  a8/5  dB/km  . (26) 

For  a frequency  of  90  Hz  with  r.m.s.  surface  roughness  of 
3 meters,  Bj^  - .0462  dB/km.  Ref.  1,  Fig.  6 (curve  3)  indicates 
that  the  attenuation  is  approximately  .043  dB/km.  Thus  to  com- 
pare the  propagation  anomaly  of  Kryazhev»  et  al1  with  the  present 
calculation  we  need  to  evaluate: 

P.A.  - TL  + 20  log1Q (r/rx)  - 9.3373  x lo"6f 3/2o8/5r  . (27) 

B.  RESULTS 

Figure  1 describes  the  quantity  J(r,z)  (see  Eq. (18a) 
as  a function  of  depth  for  two  particular  ranges.  The  function 
J(r,z)  explicitly  indicates  the  effect  of  the  rough  surface 
(through  the  amplitude  An(r)  ).  One  concludes  from  this  figure 
that  the  rough  surface  serves  to  redistribute  the  initial 
(r  • o km)  intensity  distribution  as  one  steps  out  to  long 
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(r  - 200  km)  rang®*.  It  is  clear  from  Fig.  1 that  the  intensity 
distribution  approaches  uniformity  in  the  sense  that  it  is  no 
longer  possible  to  distinguish  the  source  position  relative  to 
othir  depths.  Figure  2 numerically  verifies  equipartition 
(Eq. (8a)  ) for  two  particular  modes,  mode  1 and  mode  51.  Compar- 
ing the  approach  to  equipartition  in  the  deep  ocean  (Atlantic  or 
Pacific)  one  notes  that  the  approach  to  equipartition  in  the  Arctic 
occurs  at  a more  rapid  rate.  The  reason  is  that  in  the  Arctic 
Ocean  all  modes  are  RSR  modes  . Figure  3 displays  transmission 
loss  as  a function  of  range  for  various  surface  roughnesses.  A 
source  with  frequency  50  Hr  (60  modes  are  calculated  to  include 
all  RSR  modes)  is  placed  267  m below  the  surface.  The  receiver  is 
at  a depth  of  260  m.  As  the  roughness  of  the  surface  increases, 
so  does  transmission  loss.  For  a surface  of  r.m.s.  height  10  meters 

there  is  an  additional  loss  at  200  km  from  the  source  of  5-6  dB 
relative  to  a surface  of  r.m.e.  height  3 meters.  Comparison  of 
our  theoretical  predictions  with  experimental  observations  shown  by 

3 

H.  W.  Kutschale  explicitly  indicate  the  need  to  include  surface 
and  bottom  loss.  If  Eq. (26)  is  employed  to  account  for  ice 
attenuation,  then  the  additional  transmission  loss  corresponding 
to  r.m.s.  height  of  3 meters  at  50  Hz  is  approximately  20  dB  at 
1000  km,  in  closer  agreement  with  experimental  observation.  The 
remaining  discrepancy  can  be  attributed  to  the  bottom  loss  not  ac- 
counted for  in  our  model. 
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Figure  4 compares  the  result  of  our  formalism  to 
the  work  done  by  Kryazhev,et  al1.  The  r.m.s.  roughness  is  3 
meters.  The  correlation  length  of  the  surface  flucuations  is 
50  meters.  The  source  depth  is  50  meters  and  the  receiver  is 
located  100  meters  below  the  surface.  The  source  frequency  is 
90  Hz.  The  solid  line  represents  the  propagation  anomoly  by 
our  model.  The  crosses  correspond  to  the  work  of  Kryazhev, 
et  al1. 
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FIGURE  CAPTIONS 


Figure  1.  Mean  intensity  at  two  specific  ranges  (0,200  km) 

as  a function  of  depth,  with  the  cylindrical  spread- 
ing factor  removed.  The  source  frequency  is  50  Hz, 
the  r.m.s.  roughness  is  10  m. , the  source  depth  is 
267  m. 

Figure  2.  Verifies  the  approach  to  equipartition  of  energy 
(Eq.8a).  Same  frequency , source  depth  and  r.m.s. 
roughness  as  for  Figure  1. 

Figure  3.  Transmission  loss  as  a function  of  range  for  3 r.m.s. 

roughnesses.  Source  depth  267  m. , receiver  depth 
260  m,  frequency  50  Hz. 

Figure  4.  Propagation  anomaly  (see  Eqs.  25  and  27).  Source 
frequency  is  90  Hz,  source  depth  is  50  m.  r.m.s. 
roughness  3 ra,  correlation  length  50  m.,  receiver 
depth  100  m. 
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FIGURE  1 
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Figure  3 
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APPENDIX  A 

DOCUMENTATION  OF  WIND 
AND  PROGRAM  LISTING 

«Vr>/te  ' J OVr)  - Vr)> 

m»»n 

GAMA 

Matrix  in  master  equations  (meter) 


P 

Matrix  whose  columns  are  the  eigenvectors  of  GAMA 

DG 

Elements  are  the  eigenvalues  of  GAMA 

XKPA 

Elements  are  the  eigenvalues  corresponding  to  .the 
normal-mode  depth  functions. 

PHIN,  PHIM 

Normal-mode  depth  functions 
ZZ 

Depth  variable,  meters 
AX 

Amplitudes,  An(r) 

POWR 

Used  to  calculate  various  quantities,  such  as, 
intensity,  transmission  loss,  equipartition,  etc. 

XR 

Range  variable  in  kilometers 


. ~ " — ~r..-  •m,  ■■  j us 1 1 ■ i, 

ssss n — ^ ‘ * ' . TZ 
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CO 

Sound  speed  at  ocean  surface , meters/sec. 

NM 

Number  of  inodes 

DELTR 

Range  increment  in  meters 
NR 

Number  of  points  in  total  range  (NR-1  is  the  number 
of  divisions  of  total  range) 

SD 

Source  depth  in  meters 
KO 

Acoustic  wavenumber,  meters  1 (Real) 

SIGMA 

r.m.s.  surface  roughness,  meters 
CL 

Surface  roughness  correlation  length,  meters 

FREQ 

Acoustic  source  frequency,  Hz 

ALPHA 

Ice  attenuation  constant  (dB/km) 
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c 

c 

c 

c 


40 


PROGRAM  AlH0(  IHPUT.OUTPUT, TAPES-INPUT.TAPEd-QUTPUT.  TAPSV  TAPE  10) 
Trtla  PROGRAM  PLOTS  LOG  OF  RATIO  OF  INTENSITY  TO  POINT 
SOURCE  INTENSITY  AT  I METER  ‘ * 

PLOTS  RANGE  TIMES  INTENSITY 
PHI  IS  TO  ME  IN  RAOIANS 
COMMON/A/FREO , SO, SIGMA 
REAL  KOSQ.KO 

COMMON  CAMA( II4,I|4),P(I|4,II4) ,00< 114) ,XKPA( 114). ON ( .1  14)  DE(  114) 

,saa?a-i3i,a??,,i:g,cf,,-ac  1 u,,ax(  ,o1  • • i**  roi  > .«<  ■ o > > 

H{L/8HR«M)  /.VTL/lOrt  TLCD6)  /,VTLI/7HA(KM)-N/ 

♦.VTL2/7H  1**6* l /,HTLJ/4dZ<M)/,VTLJ/6nP/PMAX/. 

♦HTL4/3H  N /.vn.4/8H-LOG/0G// 

namelist/vausigma.  THRESH,  CL 

>*IA  CO/I  434,/,  PI/3,  14 1 5926535/,  NM/  I I 4/ , NPREV , N NM/-2  0/ 

»i£euTR/ '•5e*4/'HR/l0,/'.S0/50«/. SETA/ 1 . 6 774E-3/ 

Hew X NU  t 

3 HPREV-NPREV-NMM-2 

C SKIP  OVER  LAST  (HAVENQ)  RECORO  Op  PREV  FREQUENCY 
IFUPREV.GT.O)  RcAU(V) 
h£AQ(9)3,H,N2,KOSQ, NMM 
IF<EOF< P ) ) 77, 7 

7 FREO-CO/  SQRTCKQSQ) 

IF(NM.GT.nmm)GO  TO  76 
C SKIP  OVER  THE  EIGENFUNCTIONS 
NZl-NZ-l  • 

KO-SORT(KOSO) 

UO  40  I- I. NMM 
REA0C9) 

REAJ<9XAKPA<I),I-|,NN) 

OO  65  t-I.HM 

XKPA C I )-XKP A ( I ) **2/ ( 2 . *<0) 

JZ-H/NZ I 
OO  50  IX. HZ 
ZZtl  >-<t-l)*OZ 

reao  val 

IF(eOF(5>)  77.31 
PRINT  VAL 

IF(PrtI.EU.90.)PrtI»A5IN<  I.) 

ALPHA— 9.33  73E-6*FREQ**< I ,5)«SIGMA«*( 1.6) 

Nam I •NM—i 

64- 1 ./( I . *tiETA»H)  —3 
00  2 N-I.NMMI 
RE- I HO  9 

IP(HPREV.U.O)  GO  TO  V 
JO  6 I-I.NPRgV 
REAU<9) 

00  10  I-I.H 
SfiAJ(9) 

reaoiv)  phin 

Nl-H-I 

OO  3 M-NI.NM 
XI  ■XKPA<N  )-KKPA<  :4  ) 

A»A6S(X I ) 

XX-X  **2 
REA0(9)PHIM 

SUM-0. 5*PHIN<NZ)»PHIm(NZ)*H 
SUMA -0.3-PH IN (HZ) *PHIM(NZ) *64 
JO  60  I >2 .HZ  I 
AM- I ./( I .*6ETA*ZZ( I ) )**3 
SUmA-SUMA-PHIH< I)*PHIM(I)*A4 
SUM-SUM-PHI  H(  I)*PHIMCI  )*ZZU) 

ZMN-SUM*OZ 


63 


50 

4 


51 


10 


'•  »r- »‘¥n  " 


4MMI-SUMA*0Z 

C *RITE(6. IOOO)N.M.ZMN 

Cl  OOO  FORMAT! IX. *ZMN(*.  13 , *,*.  13.  *>  — , E I I . 4 ) 

Z2-ZMN—2 
R-XX-CL— 2 

IF(R.LE.I.c-8)G0  TO  22 
IF(R.GE.350.  )00  fo  23 
00-EXP! -R) 

00  TO  12 

22  00-1. -R 
30  TO' 12 

23  00-0. 

1 2 GAM AM , M )«CL*SIGMA— 2*DO*Z2*XX**2*KC)SQPSORT (PI) 

OAMA  ( N , M ) M 2 . *P  I ) -.CAM  AM , M) * ( I BETA-INN  I / ( XX*ZMN  ) ) ) **2 

GAMA!M.N)«GAMA!N.M) 

3 CONTINUE 

2 CONTINUE 

00  37  t-I.NM 
GSUM-0. 

OO  38  J-I.NM 
IFIJ.EO.  I )G0  TO  38 
GSUM— G5UM*GAMA( J. I ) 

58  CONTINUE 
57  GAMA!  I,  I > — GSUM 

C OO  BO  I-I.NM 

C 80  .1RITE(6.  I060X0AMA(I.J).J-I.I  > 

C *——— CALCULATES  COHERENT  DECAY  LENGTHS—— 

CALL  PLOTSoL! 275. 9HA . BEILIS) 

OO  13  I-I.NM 

13  0E(  I )— ALOO  IO(AdS(GAMA(  I,  I ))) 

00  16  J-I.NM 

16  DNU ) -FLOAT <J> 

EMAX-OEd  ) 

00  115  I-I.NM 

113  I F(OE ( I ) .GE.EMAX)EMAX-OE( I ) 

EM1N-0EC I ) 

OO  116  J-I.NM 

116  IF<OE(J).L£.EMIN)EMIN-OE(J) 

NH4-3 

NV4-8 

Xl-Fl 

CALL  GFX  (ON  ,0E.  NM.ONII  ),DN(NM).4.,6.,XI  . 

*0E(  NM ) .HTL4  ,NH4 , YTI.4 . NV4 ) 

OO  17  I-I.NM 

17  0E(I)-O. 

CALL  FRAME 

C _ 

G MAX-0. 

NMMI -NM-I 

00  35  I- 1, NMMI 

JMIN— 1*1 

OO  a6  J-JMIN.NM 

IEmP-ABS ( GAMA ( J , I )) 

IF( TeMF.GT. GMAX ) 3 MAX -TEMP 
56  CONTINUE 

53  CONTINUE 

•<R1  TE(6.600)GMAX 
600  FORMAT! *0  GMAX— .El  1.4) 

C MUST  PRINT  GAMA  ARRAY  NO*  SINCE  MAY  3E  DESTROYED  BY  TRE02 
C *RITE(6, 1062) 

C 1062  FORMA T ( *OTHE  LOwER  TRIANGLE  OF  THE  GAMMA  MATRIX  FOU.OMS*) 
OO  371  I-I.NM 
00  581  J-I.NM 

I F( A BS( GAMA (J , I > ) . LT. THRESH* JMAX )GAMA(J,I>-0. 

581  CONTINUE 

571  CONTINUE 


I •.  % j*  •«  •:?*.  rf'y: 

' HfW:.  •*, 


.1. 


CALL  TREJ2tNM,MM,0AMA,0G.fc.P) 

CALL  IMT0L2(NM,NM.DG,E,P,I£RR) 

IFt  l ERfl.NE.  0)  GO  it)  SOU 
C aRITE(6.  I059)FR£0 

Cl  05V  FORMAT(*0Irl£SE  ARE  THE  EIGENVECTORS  CORRESPONDING  TO  ACOUSTIC 
C I FREOUENCY  *,F6 .2) 

C 00  10  l-I.NM 

C70  ARl TEt  6, 1060  )(P(J,I),J-I,NM) 

1060  F0RMAT(*0*. IOEI 1 .4) 

ARI TE (6, 2223 ) FREJ.NM,  SD 

2223  FORMATt*  FREU-*,F6.2.*  * NODES-*,  1 4.*  SOEPTH-*.F7.3) 
aRITS(6, 1061) 

1061  F()RMAT(*OTHE  EIGENVALUES  ARE  *> 

ARI TE<6, 1060) (DG( I ),I-I , NN) 

JO  I a I-l.NM 

DEt  I ) — ALOGI  JtAdStDGt  I )>) 

Id  IFt I .EU.NM  > J£ ( I )-0E ( I- I ) 

Nrl4— 3 
MV4-d 
XI-Fl 

CALL  GFXtON ,0E, NM.DNt I ) .ONtNM ) .4. .6. , XI , 

*JE ( MM) . HTL4 , NH4 , VTL4,  NV4 ) 

CALL  FRAME 
JO  101  l-I.MM 

00  101  K-I.NM 
SUMd-O. 

JO  103  J-I.NM 
A I -OGt J ) *OELTR 

1 Ft  AdS(X  I ) .GE.675  .d4)G0  TO  103 
YI-P(I.J)*PC<.J)*EXP(XI) 

5UMd-SUMd*Y I 

103  CONTINUE 

C dX  IS  REPLACED  BY  GAMA  TO  MIN . SPACE  USED. 

101  GAMAt t.KI-SUMd 

C aRIIE(6. 1060) (GAMAt I. K). K-I.NM) 

N XX- l N f ( SO/DX ) * I 
a-S0/DX-FLOAT(NZA  >♦ I . 

REaIND  9 
REAJI9) 

5UMC-0. 

DO  20  J-I.NM 
REAJ(9)PHIN 

PHINT-A*PHIN(N2X*I )♦( I .-a)*PMIN(NZX) 

AXt I ,J)-PHlMf*»2 
20  SUMC-SJMC*PHI  MT*»2 

CONST- 1 ./SUMC 
aRI TE(6, 3333 ) CONST 
3333  FORMATt*  Ct)NST-*, El  1.4) 

C aRIT£(6, 3334) AX (I, I). AXt I. 50), AXt I, 100) 

JO  30  I- 1 .AM 

30  AXt I , I ) —AX ( I , I ) -CONST 

C aRI TE(6, 3334 ) AX ( I , I ),AX( 1 .50) , AXt I , I JO) 

C3334  FORMATt*  AAt I , I .50, I00)-*,3EI I .4) 

JO  d2  L-2.NR 
LMI-L-I 
JO  dl  I- I ,NM 
SUMA-O. 

DO  S3  K-I.NM 

d3  SUMA-SUMA*GAMA( I ,K)*AX(LMI ,K) 

AA(L,  D-SUMA 

C (F(MOO(L.IO).EO. I > aRITE( 6, 2222) LM I, I,  AX (L, I) 

C222  FORMATt*  AX  (*,  14,*,*.  1 4,  *)•*,  El  I .4) 

SI  CONTINUE 

82  CONTINUE 

C DO  701  L-2.NR 

C ENCJN-O. 


I 
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c 00/11  l-I.MAI 
C 7 1 I ENCl)M-ENCON-AX(L.l) 

C70I  aRITE(6, 9999  )EMCOM,L 

C 9.999  FORMAT!*  ENCOM-*,Et  1.4,*  RAMGE-40TI MES* . 1 4 . -KM* ) 

00  605  J-I.NR 

605  Xfl(  J )-FLOAT( J-l >*OELTR/l .E*3 
RH-8 

REaIMO  9 
REAJ(V) 

DO  21  M-l  .MM 

27  REA0(9)(PHI«d<N,n,I-l,MZ) 

XA-2.*AL0GI0( 12./I3.) 

XBI-AL0GI0(2.*PI/<K0*C0NSI)> 

UO  503  M-1,6 

1300"  I O0.*FLOAT<  4- 1 )♦  1 00. 

MZY-INT(QOO/OZ)*I 
A I -GOO/DZ-P  LO  A T < MZ Y ) * I . 

SUMP  1-0. 

00  31 7 M-l. MM 

PHIMA-AI*PH1MB!M,NZY*I  )♦<  I I >*PHIMB  <N,NZY) 
d I 7 SUMPI-3UMPI*AX<  I ,M)*PHIMA**2 
POARC  I > -SUMP  I 
C ARIic<6.503)POAR(l) 

C503  FORMAT!*  PO*ER(  I ,Z)-*.EI  l.*> 

POAfl!  I )— 1 0.  *! ALOG 10 ( POAR ( I >*XA*XBI  ) 

*Rl TE!6. 98 ) HOAR < I ) 

98  FORMAT!  I X,»LOC<  t/S)  ( I )**,EI  1.4) 

IF! AdS! POaR ! I ))  .LT.dO.)Pl>AR!  I )— 80. 

I F ( AdS! POaR < I ) ) .GT.  1 30.  ) POARU  ).—  1 30 . 

00  63  1*2. MR 
SUMP-0. 

OO  307  M-l. MM 

PHIMX-AI *PHIM8!  M.MZY-I )♦! I .-A  I )*PH1 MB  >M ,NZY ) 

307  SUMP-SUMP-AX!  I.M)*PrllNX**2 

POaR(  I) -SUMP 

POAR<I)-I0.  *!ALOG  10!  SUMP/I  XR!  I)*l  .6*3  ) >*XA*XB  I ) 

C POaR( t )«POaR( I ) -aLPHA*XR! I ) 

C POaR<  D-POaR!  n*20.*AL0GI01XR!i>»I.E*3*J3./l2.) 

ARlTE!6,99>l.P0ARU) 

99  F0RMAT!IX.«LGGU/S)!*,I«.*)-*.EII.4) 

IF!AdS!P!)AR!  I )).0r.30.  )P0.|R(  I )— 30. 

I F! AdSiPUAp ( I ))  .GT.  I 30.  >P0AR<  I>  — 130. 

0 IF!AdS!POAR! I >).JT. 15. )PO«R< I )-l5. 

C <*Ri  TEI  6.2323)  1 .000,  SUMP,  POaR!  1 ) 

C2323  FORMAT!*  P!*,  12  ,*.*,F7.2 .*)-* ,E  II  .4.*  LOG!  P)-*.E I I . * ) 

63  CONTINUE 

C f MAX-0. 

C 00  d27  I-l.tfR 

C82 7 IF! POAR!  { ) • GT . YMAX)  YMAX-PO aR(  I) 

C YMH-YMAX 

C 00  337  I-I.NR 

C3J7  IF!P!)MR!I).LT.YMIN)YMIN-POAR!  I) 

IF!M.GT. I )Mri--d 
XQ-0.3 

XMO— 0. 

CALL  YFXUR.POmR.NR.XHU  >,XR!  NR)  .-I  30..-30.  .HTL.MH.VfL.  IO.M.XO. 
♦XMO) 

501  CONTINUE 

CAU.  FRAME 
MH-d 

OO  601  M-1,6 

000-1 0O.*FLOAT!M“l  )*I00. 

NZY-I NT! DOO/OZ) *1 
A I -OOD/DZ ♦ I . -FLOA  r<MZY ) 

SUMP  1-0. 

OO  917  N-I.NM 


I 


!• 
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MM l NAm.n  *HM  IM6(N.  N/  T * I ) ♦ ( I • -rt  I )*HHJ  fU  IN . NZY  ) 
y I 7 aUMMI^oUMHI  *AX(  I ,M)  *KMliiA*»2 
KO..HI  I >“30491*1  .c*2 
C nRlTE<6.60J)J<)J.R')„R(  I ) 

C603  FOR«*AT(*  K<  I ,*,F/.2,*)“*.el  1.4) 

DO  73  1-2, NR 
SUMM-O. 

DO  VO 7 N-I.N4 

PHINX“Ml*PrtlNB<N.NZY*l  >♦<  I .-.1 1 >*HHIN3  CN.NZY) 

907  SUMK*SUMK*AX( 1 . 7 >*KHINX**2 
M0“R( I )“3U«M» I . E*2 
C rtRITECd, 43«3)1. D:)D.3o4P 

C4343  FORMAT! • PC *, 13 . • F 7.2,*> •*. E II . 4) 

73  CONT l NOE 
C YMa-0. 

C DO  d47  I-I.NR 

C847  IFCPOnRC  1 > .Of  .YMX  ) YMX“PO.»R(  I ) 

C nfll  l£  (6. 4 000  > Y‘4X 

C4000  FORMAT ( * M4AA**,cl  1.4) 

C DO  657  l • I , .4 R 

C857  P0“RU)“K0aR(  IJ/YMA 
IFCA.GT.  I >.4H— 6 
AO-6.5 
XMO-O. 

CALL  YFXC  XH.90.iR,  NR,  XRC I ),XS(NR)  ,0. . I . . HTL.  MH.  VTL2, 7.M.  XQ. X'AO) 
601  CONTINUE 

CALL  FRAME 
XMO-I . 

000-300. 

MZY-INTCOOO/OZ)*I . 

NX-UOD/DZ*l .-FLOATCMZY) 

OO  791  N-l.M4.l6 
NH-d 

C PHIMX-MX*PMIMBCM.MZY*I  )♦<  I .-.«X)*PHINBCN.NZY) 

OO  792  1*1 , MR 

POMri C I )*AXC I ,M)  *FLOAT(  N4) 

IFCPOrtHC  [).  or.  2.  >•<.<  I TEC  6. 26)  I , N.POMRC  I) 

26  F0R4ATC  IX,«P0nR<*,I3.*.«*.  13, •)■*.£ 1 1 .4) 

792  I FCPORRC I ) .Gi  .2 . > POMRC I ) *2. 

CALL  YFXUR.FOaH.NR.XRC I ) ,XR< NR) .0. ,2 . .HTL.NH. VTLI , 7 .N.XQ. XMO > 
791  CALL  FRAME 

OO  793  1*1. NR. 10 
XI “FLOAT (I -I  )*0ELTR/I.E*3 
IX-XI  • 

1F( I .EO.MR)CO  TO  666 

IFUX.0T.460.AN0.  C.4O0C IX, 300)  .ME.OIIGO  TCI  7V3 
666  OO  794  J-t.MZ 
3UMX-0. 

OO  795  N“I,NM 

795  SUMX“SU.4X*AX  ( I ,N  >*MHINB<  N,  J ) **2 
794  PHIMCJ)-SU4X 
YMIA-O. 

00  89  J“I,NZ 

89  IF(PHIN(J).GE.YMIX)YMIX-PHIN<J) 

00  91  J“l ,NZ 
91  PHIN(J>-PHlN< J)“l .E*2 
Nri3*4 
NV3-6 

CALL  GFX ( ZZ . PHI M , NZ . ZZ ( I ). ZZC  NZ > ,0. , I . , XI . YMIX. 
♦HTL3,Mrt3,vrL2.NV3) 

CALL  FRAME 
793  CONTINUE 

CALL  PLOTCO. 0.0. 0.999) 

OO  TO  4 

500  “HI  TEC  6,  1 06  J)  IE  HR 
1050  FOH4AT(*OIERH-  *» 120) 
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f 


I 


ierr-ierr-i 

I FI IERR. GE.  I.  AND.IERR.LT.  25 I )«RITEl  6,  1060)  (DC(  I ) . I- 1 , 1ERR) 
GO  TO  5 

76  -RITE16,  I080JNMM 

1080  FORMATl*OERROR.()NLY*,  I5,-*MODES  ON  DISK*) 

77  SWF 
ENO 

C 

SUBROUTINE  IMTOUINM.N.O.E.Z.  IERR) 

INI  EGER  l .J  •<  >L  «M.N , 1 1 ,NM, MML  .IERR 
REAL  0<N),E<N),ZINM,N) 

REAL  8,C,F,0,P,R, S, MAChEP 
MACHEP  - 2. **(-47) 

C 

IERR  • 0 

IF  IN  .EQ.  I)  GO  TO  1001 
C 

00  1 00  I - 2.  N 
100  E ( I-l  ) - Ed) 

C 

E(N ) - 0.0 
C 

00  240  L - I.  N 
J ■ 0 

C «««««»««  LOOK  FOR  SMALL  SUd-OlAGONAL  ELEMENT  ********** 

I OS  JO  110  it  * L,  N 

IF  IM  .EG.  N>  GO  TO  I 20 

IF  IA8SIEIM))  ,LE.  .MACHEP  * lAdSIDIM))  * AdSlOl M*l ) ) ) ) 
A GO  W 120 

110  CONTINUE 

c 

120  P - OIL) 

IF  (M  .EQ.  L)  GO  TO  240 
IF  (J  .EO.  30)  GO  TO  1000 
J - J ♦ I 

C **********  FORM  SHIFT  ********** 

G - CO  CL*  f ) - P)  / <2.0  * £<L7) 
k - SQRi  IG*G*1  .0) 

0 • JIM)  - P ♦ £<L)  / (G  * SIGN(R.G)) 

S - 1.0 
C - 1 .0 
F - J.O 
MML  • M - L 

C **********  for  [««*- 1 5TEF  -I  UNTIL  L DO  — ********** 

JO  200  II  • I . MML 
I - M - II 
F • 5 * Ell) 

S • C * El  I ) 

IF  UdStF)  .LT.  AdSl G) ) GO  TO  150 
C ■ G / F 
R - SGRTtC*C* I .0) 

£11*1 ) • F * 2 
3 - 1.0/  R 

c • c * s 

30  TO  160 
150  5 - F / G 

R • SORTtS*S*l .0) 

et  i*i  > • g * r 

C - 1.0  / R 

s - s * c 

160  G ■ 011*1  ) “ P 

R - toil)  - G)  * S * 2.0  * C * 8 
F • S * R 
3(1*1)  • 0 ♦ P 
G • C * R - a 

c **********  FORM  VECTOR  ********** 


i 
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t 

* 


130 

200 


JO  130  K - I.  N 
F • Z(K, 1*1 ) 

2( K,  I*  I ) - S * Z( K,  I ) ♦ C * F 
Z(K.I)  ■ C * Z(K.I)  - S * F 
CONTI NJE 

CONTINUE 


3(L>  - 0<L)  - P 
E(L)  - J 
c(4)  ■ 0.0 

jo  ro  105 

2*0  CONTINUE 

.*»»•»*•*»  OH0ER  EIGENVALUES  AND  EIGENVECTORS 
DO  300  II  - 2.  N 

[ * II  - 1 

.<  - I 

P ■ D(I) 


•JO  260  J - It.  N 

IF  (D(J)  .GE.  P)  GO  TO  260 
K • J 
P - OCJ) 

CONTINUE 


260 


230 


IF  (K  .£0.  I)  GO  TO  300 
0< X)  *3(1) 

3(1)  * R 

30  230  J • 1 , N 
P - Z(J,I  ) 

Z(J,I>  - Z ( J . < ) 

ZU.K ) - P 
CONTINUE 


300  CONTINUE 


JO  TO  1001 

C *•**»•«*»*  SET  ERROR  — NO  CONVERGENCE  TO  AN 
C EIGENVALUE  AFTER  30  ITERATIONS 

1030  I ESN  * L 

•iSl  i'e(  5.2030)  L 
2360  l-OW4Al(*OL*  *,123) 

1 JO  I REiURN 

C *•*»*•**•*  LASi  CARO  OF  IMTQL2  ********** 

ENJ 

SUBROUTINE  TREC2I  NM,N,A,J,£,Z) 

W 

integer  i.j.k.l.n,ii,nm,jpi 

REAL  A(NM.N) ,D(N) ,E ( N ) , 2 ( NM, N ) 

REAL  F.G.H.rtrl. SCALE 
C REAL  SORT, Ads, SIGN 

C THIS  SJdROUTINE  IS  A TRANSLATION  OF  THE  ALGOL  PROCEDURE  TRED2, 

C NUN.  MATH.  II,  131-195(1963)  BY  MARTIN,  REINSCH,  AND  WILKINSON. 

C HANJdOOK  FOR  AUTO.  COMP.,  VOL.  I I-LINEAR  ALGEBRA.  2 1 2— 226( 1971 ) . 

THIS  SUBROUTINE  SEDUCES  A SEAL  SYMMETRIC  MATRIX  TO  A 
SYMMETRIC  iRIDIAGONAL  MATRIX  USING  ANO  ACCUMULATING 
ORTrlOGONAL  SIMILARITY  TRANSFORMA TIONS. 

ON  INPUT- 

<M  MUST  dS  SET  TO  THE  ROn  DIMENSION  OF  TWO-DIMENSIONAL 
ARRAY  PARAMETERS  AS  DECLARED  IN  THE  CALLING  PROGRAM 


C 

C 

C 


c 

c 

c 


. 

■ » II  ■ - ■ 1-  « — 
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DIMENSION  STATEMENT, 

N 15  fHE  ORDER  OF  THE  MATRIX. 

A CONTAINS  THE  SEAL  SYMMETRIC  INPUT  MATRIX.  ONLY  THE 
LO.iER  TRIANGLE  OF  THE  MATRIX  NEED  EE  SUPPLIED. 

ON  OUTPUT- 

J CONTAINS  THE  DIAGONAL  ELEMENTS  OF  THE  TH IDI AGONAL  MATRIX. 

E CONTAINS  THE  SUdDUGGNAL  ELEMENTS  OF  THE  TSIDIAGONAL 
MATRIX  IN  ITS  LAST  N-l  POSITIONS.  Ed)  IS  SET  TO  ZERO. 

Z CONTAINS  THE  ORTHOGONAL  TRANSFORMATION  MATRIX 
PRODUCED  IN  THE  SEDUCTION, 

A AND  Z MAY  COINCIDE.  IF  OISTINCT.  A IS  UNALTERED. 

QUESTIONS  AND  COMMENTS  SHOULD  SE  DIRECTED  TO  8.  S.  GASdO«. 
APPLIED  MATHEMATICS  DIVISION,  ARCONNc  NATIONAL  LABORATORY 


DO  IOO  I » I.  N 

JO  1 00  J • I , I 
ZU.J)  • A < I . J > 

CONTINUE 

IF  ( I .EQ.  I)  GO  TO  320 

**»*»*♦»**  PUR  I»N  STEP  -I  UNTIL  2 DO  — ********** 

JO  300  II  - 2,  H 
I * N ♦ 2 - II 
L ■ I - I 
n ■ j.o 
scale  *0.0 

IF  (L  .LT.  2)  GO  TO  130 

•»•*»****»  SCALE  RON  (ALGOL  TOL  THEN  NOT  NEEDED)  ********** 
JO  1 20  < » I , L 
SCALE  * SCALE  * A dS ( Z (I , X ) ) 

IF  (SCALE  .NE.  0.0 ) GO  TO  140 
id)  - Z(I.L) 

JO  TO  290 

JO  IdO  < • I , L 

Zd.X)  - Zd.O  / SCALE 
H * H * Zd.X)  * Zd.X) 

CONTINUE 

r - Z ( I ,L) 

J - -SIGN!  SDR  T(H)  ,F) 
s(I)  ■ SCALE  * G 
H - H - F * 3 
KI.L)  • r - 3 

■ J.O 

J.)  240  J « I , L 

Z(J.I  ) - Zd.J)  / (SCALE  * H) 
j m 0 . D 

**********  FORM  ELEMENT  OF  A*U  ********** 

DO  180  X - I,  J 
0*0*  Z(J,X)  • Zd.X) 

JPI  ■ J ♦ 1 
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IF  (L 

.LT.  JPI)  GO  TO  220 

DO  200 

< ■ JPI.  L 

200 

G ■ G 

♦ Z«,J)  * ZU.K) 

**********  F 

ORM  ELEMENT  OF  P ********** 

220 

HU)  - 

G / H 

F - F 

♦ EIJ)  * ZCI.J) 

240 

CONTINUE 

HH  * F / (H  ♦ H) 

**********  FORM  REDUCED  A ********** 

30  260  J - I,  L 
r - ZCI.J) 

D*EU)-HH»F 
HU)  * 3 


i 


C 

C 


c 


c 


DO  260  X * 1.  J 

2CJ..O  - ZCJ..O  - F * ECK)  - G * Z(I,K) 

2 60  CONTINUE 

JO  260  K » I . L 

230  Z<I.K)  » 3 CALS  * ZCI..O 

290  0<  I)  * H 

300  CONTINJS 

320  0(13  * 0.0 

E(  I ) - 0.0 

**********  ACCUMULATION  OF  TRANSFORMATION  MATRICES  ********** 
jo  600  r » I.  .>1 
l * I - I 

IF  <0<I)  .HO.  O.J)  GO  TO  330 

JO  360  J - I , L 
J * 0.0 

DO  3*0  < * I . L 

340  j « G ♦ ZU.O  * Z(K.J) 

JO  360  K • I.  L 

Z(K.J)  * IIK.J)  - G * ZUC.I) 

360  CONTI NUE 


360  J(  I)  • 2U.II 

Z(  I.I)  - 1.0 
IF  ( L .LT.  I ) GO  TO  500 

JO  400  J - I . L 
ZCI.J)  * 0.0 
ZU.I  ) » 0.0 
400  CONTINUE 


300  CONTINUE 

£ 

RETURN 

C **********  LAST  CARO  OF  TRED2  ********** 

END 

SUBROUTINE  7FX(  X , f .N.XMIN.XMAX, YMIN, YMAX  JNTL.NH,  VTL.NV  ,M, XQ.XMQ) 
J I MENS  ION  TEXT!  7)  ,X(,I)  , Y(N) 

COmNON/A/FREJ. SO.  SIG.4A 

DATA  (TEXT!  I ),!■!  ,7)/IOrtOO*  I 00. 00.  I OHOO*  200.00, 

♦ IOHOD*  3 00 . 00 , 1 OHOO*  400.00, 1 OHOO*  500.00, I0H0D*  600.00, 

♦ I OHOO-  700.00/ 

•<R[  t £(6,201  4)  (X(  I ),  1*1  ,N) 

•iRI  l£( 6, 20 1 4 ) (Y  ( J ) , J*  I ,.N  ) 

2014  FOR  *AT( IX, IOHI I .4) 
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l 


IF<NH.L=.0)G0  it!  10 
OX«<  XMAX-xMIN  )/5. 

JY-C  Y.MAX-YMlNJ/5. 

XB»0.3 

ifcxo.ej.o.3) xa^a.s 
I P<  <a.£Q.0.5)NH— NH 

CALL  AXI5(0.5,Xa.Hi'L.NH,5.  ,0.  .XMIN.OX  ) 

CALL  AXIS(0.5.0.5.VrL.NV,3.,90.,YMIN,DY) 

CALL  SYMBOL ( 0.5 ,XQ, . 1 ,3H  p-,0.0,3) 

CALL  NUM6ESC0.6.X0,.  I ,FR EO.O.O,  I ) 

CALL  SYMBOL<2.J.AO, . I .3HSU-.0 .0, 3 ) 

CALL  NUMdER<2.3.XO..  I ,50.0.0, I ) 

CALL  SYMBOL (3 .3 ,XU, . I .6HSIGMA-.0. 0.3) 

CALL  NJMdSR<3.rf.A0. . I .SIGMA, 0.0, I ) 

I rCXMO.EQ.  I • )2»M 

I F(  AMO. EO.  I . ) CALL  5YM30L(  3 .6,  6.0,  . I , srIMOUE  .0.0,3) 
I r ( a MO . EO . I . 1CALL  NUMBER ( 4.0, 6.0, . i ,2,0.0, I ) 

10  CALL  LIMEC A , Y ,N , I .O.O.XMl.'l— OA/2. , DX , YMI M— OY/2 . ,0Y) 
CETA-C  Y(  J)-YMIN)/0Y*.3 
lr( Y ( N) .LT.O. >2E TA-aB-ABSC Y< NJ-YMAX )/DY 
C XciA  OEFENOS  0.M  CHOICE  OF  AXIS 

C RE-AJJUSf  AS  SEOUIREO 

IF( XMQ.EO . 1 . )G0  ro  12 

CALL  SYMBOL (X )/DX».6,Lci A, . OB.TEXTC  M) .0., 10) 

12  RETURN 
Ha  0 

SUBROUTINE  GFXC  X , Y.N.XMI N, XMAX, YMIN , YMAX, X I , YMI X . 
*ri  7L  3 , NH3 , V Y L3 ,MV  3 > 

0 1 MENS  ION  X ( N ) . Y ( N ) 

COmmCWA/FR EO , 50 , S I GMA 
.••SI  TECS,  401  A)  C A ( I ) , I » I , N) 

. NIXSC6, 4014) CYC J),J»I ,N) 

4014  FOR  aaVC  lX.lOel  1.4) 

IFCXI.  )E.FI  ).<SirE(6,4024)Al,  YMIX 
4024  FOR  4Al"(  I X,*FMAx(*,F6.  I ,»  <«)•*, Ell. 4) 

IFUH3.LE.0) JO  TO  10 
)X»<  X MAX— AM  IN  )/5. 

JY«t  YmaX-Ym I N )/5. 

TALL  AX  1 5(0 .3 ,0.3,riTL3,-iNri3, 3 . , 0.  , AM  IN, OX) 

CALL  AAIS<0.3,3.3,vrL3,fW3,3.  .90. . YMIN.DY) 

CALL  SYMBOL <0.3, 5. 3,.  I , 3ri  F-.0.0.3) 

TALL  )U.M5ER(0.d,S.3, . I .rREO.O.O,  1 ) 

TALL  SYMBOL <2. 0,6. 3,. I , 3HS0», 0.0, 3) 

CALL  NUMBER  (2  .3 , 6 ,3, . I ,50.0.0,  I ) 

CALL  SYMBOL ( 3 .3 , 6.5, . I .6HSI0MA-, 0.0,3 ) 

CALL  NUMBER (3. 8, 6. 3,. 1, 5 1 GMA, 0. 0, I) 

IFC  xI.NE.FI >CALL  SYMBOLC  2.0, 6.0, . 1 . 6HRC  KM ) » , 0.0, 6) 
IFCAI.N6.FI >CALL  SU  ABERC  2. 3, 5 .0 , . l.AI.O.O, I ) 

10  CALL  LINeC  X « Y , I ,0,0,AMI.<— OX/2  • , OX,  Y Ml N-OY/2 . , OY) 
SE TJ  SIM 
e.»J 
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APPENDIX  B 

DOCUMENTATION  OF  NORMAL  MODE 
PROGRAM  FOR  THE  BILINEAR  PROFILE 
AND  PROGRAM  LISTING 
P 

Normal  modes  $n  as  a function  of  depth 
U 

Eigenvalues  un 

STK  

Related  to  eigenvalues  via  STK  =*  kQ  * l-un 

where  kQ  is  the  source  wavenumber 

2M  Depth  coordinate  (variable) 

ZMAT 

Depth  coordinate  at  which  up  and  down  integrations  meet 

ZMEET 

Depth  coordinate  where  the  two  velocity  gradients  meet 
EN 

Index  of  refraction,  n(z),  as  a function  of  depth. 

VM 

The  potential,  l-n2(z),  in  the  normal -mode  equation 

DELZ 

Depth  increment 
KO 

Source  wavenumber 

KOSQ 

Source  wavenumber  squared 
BC 

Boundary  condition  value  for  eigenfunction 

SLOPE 

Slope  of  eigenfunction  at  the  boundary 

TPT 

Turning  point  where  eigenfunctions  turn  from  oscillatory 
(exponentially  decay)  to  exponentially  decaying  (oscillatory) 

H 


Ocean  depth  (m) 
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FBGN 

Source  frequency  (HZ) 

IM 

First  mode  calculated  (also  equals  the  number  of  zero  cross- 
ings in  the  mode) 

NM 

Last  mode  calculated  (equals  number;,  pf  zero  crossings  in  the  mode) 

NMINC 

Increment  of  successive  modes  to  be  calculated  between  IM 
and  NM. 

NZW 

Number  of  depth  points  (spaced  at  DELZ)  at  which  the  normal 
mode  depth  functions  are  calculated  and  stored 
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USER  .'(UMBER i adit 
BEGIN  TEXT  EDITING. 

7 ll* 

PROGRAM  MAINdNPUT.OUTPUT.  TAP  E5- INPUT,  TAPE6-OUTPUT.TAPE9 ) 

C INSTRUCTIONS  FOR  CHANGING  DIMENSIONS 

C CHANGE  THE  DIM  IN  U AND  SIX  TO  # OF  MOOES  REG 
C STEP  I.  SELECT  NEn  DIM  FOR  P,  -IS  TO  20  PCT  ABOVE  ZM,  VM 
C STEP  2,  as  /OLD  OIM/./NEn  0IM/I2 

C STEP  3.  SELECT  NEN  DIM  FOR  ZM.  VM  - 2**N*I  FOR  SOME  N 
C STEP  «.  RS  /OLD/./NEa/I 10 
- C STEP  5.  RS  /OLD»P/./NEn*^/ 

' DIMENSION  U(20) ,P(960I ).STK(20) 

COMMON/GEN/UG .2MI 800 1 ) .ZMAT. F ACT. DELZ . XO .01 VCX. SQ8. SUB.  KOSO.  VM( 
♦8QQi ) .H12. Hid. EXPMAA.OSCMAX.dC.TPT, SLOPE. EPS, B 
REAL  JC0.K030 

DATA  CO, ZO, B. G/ 1 434 . , 0. , 1 300. ,.017/ , EXPCT. OSPCT/ . I , . I / 

OATA  ZMAT/I05./.ZMEET/337.  3323*27/ 

OATA  NPTS.PI .QIVCX.UTEST/  800 1 .3. 14159265358979. I .£-40. 1 .E-4/ 
DATA  H.nX.PCTERR.NP)O/4000.,0.5.0.0I.960I/.SLINIT/I.E-I0Q/ 
f NAMELIST/VAL/FBGN.FMAG.NFQ. IM.NM.NMINC, NZN 

C EN(Z>-I ./( I .♦EPS*<2.*<Z-ZD)/8“l.*EXP(2.*(ZO-Z)/8  ) )) 

EN(£>-l./<l.006l0d525-M. 67741  E-5*Z) 

V<Z)-I.-£NCZ)**2 

REnINO  9 ........ 

EPS-3C/CO0.3 

, EPSQ-SQRT CEPS ) . . 

DELZ-K/  ( NPTS—t ) 

H12-0ELZ**2/I2.  - - 

, HI6-I6./21./0ELZ 

SQ8-.9  . 

FACT-PCTERR4240./H 

DO  10  l-I.NPTS  - - - . 

ZMU)-OELZ*U-t ) 

IF(ZMtl)  -GT.ZMEEDGO  TO  300-  - — . . 

VM< I > -I .-< I ./ Ct . *-3  «48765E-5*ZM(  l > ) ) **2 

GO  CO  10  - 

300  VMLI>-V<ZMttJ) 

ia  CONTINUE  — : -Z.  . 

IMlO-FLOATtNPTSI-ZMAT/H  . 

IFtZJtt  IMIQ-H ) -ZMAT.LT. ZMAT-D4C IMID7V  IMID-IMID*!  

NMIO-NPTS-IMID*l  - •-  / 

IMIOO-IMIO  — . - ..  

NMIOO-NMIO 

l FLAG-0  --  . , ' 

OOi-IMIO-l 

NMOI-NMID-I  - - 

XSHFT -NPMX-NPTS 

vlept-vmi  i) 

VRT-V (H) 

ARITEt6.8DVLEFT.VRT 

81  PORMATCIX.*VL£FT-*.FI0.5.*  VRT—.FIOJ) 

KOSO-I. 

1 READ  val 

IF<EOF(5>> 100.2  ^ 

2 PRINT  VAL 
IFtFBON.LT.O.J  STOP 
NSXP-< NPT5- I ) /( NZN- D 


IM-IM4I  - . 

NtHNM*l 

FRAT-IO.**< l./FMAG) 
IFtFBGN.EQ.O. ) GO  TO  23 
r-FBGN/FRAT 
23  DO  22  II-I.NFO 


( 
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n«2.<»Pl«F 

. OUKSOKOSQ  “ 

ko-a/co 

K0SG«K0*K0  "■ 

ARITEC9)  B.H.MZrt.KOSO.NM 
ARITSC6.3545) 

3545  FORMAT l *nROTE  HEA06R*) 

RATIO-K030/OL0KSQ 
00  II  K-I.MPTS 
II  VM!K)»RATIO*VM(IC) 

A-B8-ZMAT 

NRR-0 

KNTM-0 

00  20  I-IM.NM.MMINC 
M-I-l 

IFtIFLAG.EO.O)  00  TO  15 

IMI0-I4ID0 

NMIO-MM(DO 

IMDt-lMIOl 

NMDI-NMID-I 

IFUG-O  -•-= 

15  AK8»lAK*M)*Pl 
IFCKNTM-I ) 16.17.18 

16  UG»!C3.*Pt*3.48765&-5/!2..*K0)  l*!FU)ATI2*M)*l  .5)>*«<2./3.) 
IFU. 487656-6. LT.OIVCK)  UG«CFL0ATC2*M-M  )*PI*.5/H/KO)**2 
OIMC-UG 

00  TO  13 

17  OIHC-U(l).*CtCFLOAT(24M)»l.5)/(FLOATC2-*IM>-.5))**t2./3.3-l. 
GO  fO  19 

IS-  UIMOUCKNTM)-U!KMTM-l>  — - --  ---  - 

19  IM*»UCIWTM)»UINC 

13  UCL-0.  - 

IF (KNTM.GT.O ) UGL-OCKNTM) 

DO  50  U.-I.I0  - • - •. . •.-"-i-f  \ - 

A— I.E50  ; . 

- 08-UE5O  — — ' — • * r 

C MRITEX6.72)  ZMAT.H.BB  - . - 

IFllJG.LT.VLEFD  CAIX  TURMCO..Z0.A)  - 

IFIUCXT.VRT)  CALL  TURK CZO.H.BB) 

C aRITEI6»72)ZMAT.H»B8  > 

C 72  FORMAT! IX.*ZMAT,H.aB»*.3FI0.3) 

iF<liC.GE-.VL£FT.OR.UG.G£.VRI)  GO  TO  60  - 

XIMT-O. 

UD6R*0.  . ...  ..  _ * — — . . 

KL-FUMTIMPTS-M  )*A/H»I . 

KU-FLOAT!  NPTS-I  )-*88/H»  l ► — - 

00  40  KX-KL.KU 

- ARG-OO-VACHO/KOSQ  • 

IF4AR0U.T.0. ) MR ITE! 6.26 ) ARG 

26  FORMAT!  I X.*ARO«».SI  1.4.)  - .„  A-Lji-  * 

IFIARG.LS.OIVCX)  GO  TO  4a  .......  . . 1 * 

SO-SORTCARG)  - — - — — • 

UOER-UOER-K3.3/SO 

XIMI-XINT*Sa  * 

40  C0KTIMU6 

XIfiC*KO*xrNT*OELZ  - *.* — — •>  '. 

UOE«-ICO*UDER*OELZ  ' • .?  • 

IPCABStXIME— <IB>.LE^)^JHWia>  GO- T0-60 
C 4RITE!6,25)U0ER 

C»  FORMAT!  IX.  «il06R-».et  1 .4)  

UINC"! AK8~X  MT)  /UDCR. 

UG-OG^JINC 

50  tPCUO.LS.UOD  UG*OGL*.  I*A£SCU  IMC)  ..  . 

aRlTBId. 1002)  M.XIMT.4K8  _ • . 

1002  FORMAT !*OMU  4KB  COHV  FOR  MODE*.  1 6,  •ACTION**.  El  I .*»*4KB**.I 


/.>  *■  r 


**  • 


lit . 
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DO  70  KX-1.50 
SUU-KOSO-UG 
IrUCK.EQ.  I ) <30  TO  51 
A— I.E30 
a a- i .£5o 

XFXUO.LT.VLEFT)  CALL  TURN<O..ZO.A>  - - - 

IFCUOXT.VRT)  CALL  TURNCZU.H.08) 

C ARITS(d.72>ZMAr.H.3S 

51  ZL-AMAXIXO..A) 

zr-amxnu88.h>  • 

OSCMAX-XZB-ZLJ/FLOATl  D-USPCT 

EXPN AX-O.  " ; - 

IF(A.M6.-I.E50>  EXPMAX-ZL-SXPCT 

52  BC-O. 

SL0PE-SLINIT*(-1 . >**M 

cam.  MUMERXP.MPMX.I  .KM.PINTI.NOOUIMXDI 
SXPMAX-O. 

rpr-d8  __ 

XFXBB.NE.I.E50)  EXPMAX-XH-ZR)*SXPCX 
56  - BC-SLXNIT 
SLOPE-O. 

CALL  tiimem P< KM* I ).NPMX-KM.NPTS.KDUM.PINT2,N002.NMID) 

IREF-KM-NMIQ 

XAB-KM-KDUM 

CALL  IMTERPCIMXO*MMIO-KDUM.KM.P.NPMX.O> 

CALL  INTERP1 IREF*IMID-KM,  (Ada  P,NPMX,KM*NPTS*I  > 

SUM-0. 

DO  63  K-IMOI.KM 

XFXP(K)— P1K-I ) .LT.O..OR. P(K) .EQ.O. > NXX3I-N00I-I 
63  SUM-SUM— P<K  > **2 

P I MX l -P I MT I -0 ELZ- ( S UM-O . 5* ( P U M I 0 > —2 -P t KM ) *»2 ) > 

XRFI-IREF* 1 . . 

SUM-0.  -•  • . ;•  . . t,  . ■ 

OO  66  K-IRFI.IAB  ....  , . . ...  V V V ' - 

SUM- SUM-PI  K)  —2  — ----  ' 

66  IF(PIK)mP<K-1 1 . Lr.O.. OR. P<K). £0.0.1  NO02-M002-1 

' N001 -N001  ♦N002  % 

P IM X2-P  1MT2-0ELZ*( SUM-O. 5* ( PU REF ) «*2-B (I AB ) —2 J ) 

IFlP-lIMID)*PlIREF>.GT.O.>  GO  TO  67  

iFUG-tFLAO*i  - 

C - XF(M.E0.0>«*RXTEC6.247)M,IFLA0  - - — — 

C 247  FORMATl*ITERS.  FOR  COKVERO.  IMOOE  *,I3.—*.I3> 
tFUFLAG.GT.SOl  CO  TO  78 
KL-OSCM  AX/D  ELZ 

IFUMDI.GT.KM1  CO  TO  246  - 

OO  240  K-IMOI .KM 

lF(P(K)*PCIREF-{MIO-K).cr.OXWOCl  CO  TO  260  - -- 

240  CONTXMUE 

24*  IFURF1  .CT.IABI  OO  TO  252  - - - - — — - - 

OO  250  K-XRFI ,IA8 

XFtP.aCl*PUMlD-IRfiF-K>.Gr.OIVaC>  GO  C0~270  - - 

250  CONTINUE  • 

254  IFXMOOI.EO.M1  CO  TO- 996  . — — - ' 

OO  TO  61  , ~ J:  - ..L 

260  XMXO-K  - - — ~V-  - 

IFXK-KL.CT.KM)  KL-KMHC  ~ - 

IFIPOt-KLI-PX  IREF*XMIOHCHCU.GE.QIVCXl  XMID-XMIO-KL 
GO  CO  28Q  ....  . A w* 

270  XMXO-XMIO-KMHEB  - . .... 

IFXK-KL.CT.IAB)  KL-lABHC  ' " 

- IFCF(K*KL1*PCIRSF— XMXO-K— KLl.GCLOIVCd  IMID-XMIOHCL 
280  MMXD-NPTS-lMXO-l  • - • ...  . ' . - ' 

UCI-XMID-I  ...  ---r--"  . | ■■  - ••••••—  -■ 

NIGI-NMIO-I  , J-  *•'» 


\ 


■ ...«••».  I 

a..  • 1 

— -T 1 


Ty\v~!  -jr-  r~r— — 

<-*>4  ...  / ./..y  u?  r ^ 1 
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• 47 

48 


e c 
c 

T cc 
c 

- C24J 


IPU00I-M1  7l^«^75  _ 

RATIO*P<lMID)/P<IREF>  • .--  r " 7J 

CALL  06ftIVtP.IMlO.IMlO. I .QERL) 

CALL  OERlVtP.IREF.lMlD.-I.DERR) 

UIi*C»PlINID>*<OERl/RAriO-OERR)/<<PtNTI/RATro*PINT2*RATIO>*KOSO) 
IFCUG*UINC.OT.UGH)  UlNC-0.5*(  UGH-JG) 

IFtUG*UINC.LT.UOL)  UINC-0.5*<UGL-UG) 

OO  TO  79 

UOL-AMAXKUOL.UG) 

OO  TO  77  •- 

UGH*AMINl(UGH,UO> 

UINC*0.5*(UGL*UGH>-U0 
IP(m.LE.O.OR.NOiJI.LE.O)  OO  TO  79 

uoer*ug*<fluatc  m>  /float < noo i >-i . > 

lP<AaS<UlNC).GT.2.*AflS<UDER)>  D INC -ODER 

UOUG*UINC  

IF<A3S(UlNC).LE.UTESr*ijG>  CO  TO  ao 
CONTINUE 

«RICE<4.II70>  M.iCK,  I FUG 

> FORMATI IX.*CONV  FAILURE  FOR*.  13.*  EVALUE*.*KX«*.  13. * IFLAG**.I3> 
STOP 

DO  82  K-I.NMDI 
PP-?(ICM*NMID-K) 

IF(PP.E0.I.E6O>  P(NPMX-K*I)-PP 
IF(PP.NE.  I.EdQ)  P(NPMX-K*I  )-PP*RATI<) 

XL- 1 MO  I 

KU-MINOC KL-KSHFT- I . NPTS) 

OO  44  K-KL.KU 
P (K) *P ( NPMX-K*!CL ) 

IF(KU.GE.NPTS)  GO  TO  90  - 

KMOV— NPTS-KU 

OO  34  r-I.KMOV  - 

P < NPMX-X*  I ) *P (NPTSHC*  I ) 

xl-xu-i  ... . _ , _ 

GO  Ta  83 

lFOJQ.OS~l.Oy  OO  TO  2»  — — 

CALL  IMTERPC I ,NPTS.P.NPTS*0) 

PINT1-0.  - . - . L ...  . . . 

00  92  K*I.NPTS,NSXP  . :.  . * 

P INTI  *P  INTI  *P(K)**2  . _ 

PINTI-SCRTC OELZ*FLOAr<NSKP)*( PINT  1-0. 5*1  PCI  )**2*P<NPTS) **2  > ) ) 

00  91  K*l  .NPTS.NSKP  . . 

P<K)*P(0/PINTI 
IFC.iZN.Lc.  I ) CD  TO  94 
*RITE(9)  (P(K).K-I. NPTS.NSKP) 

aRITEI 4.3548)  NZ4.I.A.B8  ... 

i FORMAT t *OMROTE  EFCX*.I4.I5.*  TPTS-*.2FI0.I  ) 

KNZM*KMTM*i  — _ ........ 

U(KNTM)-UC 

lp<  UO.LI.VLEFT.  AND.  UO.  LX.  VST)  NRR-NRR-I  ...  .... 

ST1CIXNTM)-X0*S0RTC  I ,-UG)  . ... 

*RITS(4*  7070)  KNTNrSTKC  JCNTM)  

FORMAT! IX.*STK(*»I4,*)«*,FIQ.7> 

CONTINUE  — ’ - 

«RITE(91  (STK(I).I-I.KNTM) 

NRITE<4»4448)  F,NRR.XNn«*(ST)C(I).I*».)a«I»i  ... 

FORMATI  F8.3  .213*  I OEM.  5)  . ...  . 

CONTINUE  - — • • - -i’A'r-. -w-  • • ...  , 

0IF23-STKC2)-8nC(3)  "•  - * •• 

0IF2S24*«TK<2»>-STK(26F  .-A 

OIF9Q9(*STIC(90)-STK(9I)  ■*.  . 

QIF202I*ST1C(  120)— 5TKC  1.21)  ' * ■ ' ' - . 

aRITE( 4.243 ) 0 IF23 .OIF2524. OIF909 I .0 IF2021 

FORMA  r<  45 1 1.31  - 

oo  to  i ■' • - ' - - l • - 
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'\ 


1991  FORMAT  <*S  I GN  CHANGE  AT  NE*  MATCH  PT,  MO0E-.14)  ' - 

GO  TO  646  ~ 

994  NRITE(4.I994>  P!!MtD),P! I REF) ,IM ID. I REF, KM. I A3 

1994  FORMAT; *PRO0LBl.  NOOE  AT  MI0*.2El I .3.4I6J 
446  aRITE(4,6648>  F.NRR.KNTM. (STXCI ) .1-1  .KnTM) 

100  STOP  , - 

ENO 

SUBROUTINE  NU*R<P.NP.IS.ICMAX,PINT.NO0ES,IMID) 

DIMENSION  P(l) 

REAL  KO.KOSO 

COMMON/GEN/ UG.ZM(  800 1 ) .20. FACT. DEU, KO. DI VCK. 308, SUB, 

♦ KOSO.VMiBOOl >, Ki2.HI4.EXPMAX.OSCMAX.PS.TPT.SLOPE.EPS. 8 
NOOES-O 


I- 


10 


l 


e 

c 

c . 

t‘ 

c 


120 


ISO 


OO  10  1*1  ,NP 
P(I)-I.E40 
PINT-0. 

ZS-ZM(IS) 

Fl-VMI IS  >-SU8 

IF(OSCMAX.LT.06U)nRITE!4, 1 61 )OSCMAX.0EU 
IAI  FORMAT!  IX, -OSCMAX—.FI0. 5.*  0ELZ—.FI0.5.*  II*) 

IF(OSCMAX.LT.DEU)  GO  TO  300 
HMAX-AMAX  I ( EXPMAX  .DEU) 

I F< HMAX .EO. OSCMAX ) OSCNAX-I .00001 *OSCMAX 
IF<A8S<TPT).EQ. I.E50)  HMAX-OSCMAX 
IF(A8S(FI l.Lfi.OlvaC)  HTST-HMAX 

IF(A8S(FI ) .GT.DIVCX)  HTST-AMI Nl (S08*EXP( ,2*AL0G( FACT/A8SC FI )**3 
* ) > .HMAX) 

IF(rirST.LT.0EU)4RITE<6.  I62)HTST,DEU 
142  FORMAT! IX, *HTST«*. FI 0.5, * DEU-*.FI0.5.*  18*) 

- IFiHTST.LT.OEU>  00  TO  300 
ZZ-QELZ. 

00  120  NSK-I  . 100  ' .... 

ZZ-2.*ZZ 

IFiHTST.LT.Z2T)  GO  TO  130  ..... 

CONTINUE 
MRITEI4.I020)  HTST 

1020  FORMAT! *MAX  STEP  SIZE-*, El 0.3. * TOO  BIG*) 

STOP  . 

HH— 2Z/2. 

N ASK-2**!  NSIC-I)  . . 

FF-HH*)H/I2. 

NSK-NASK*SIG.N(  I .O.ZO-ZS)  — 

Pil)-PS 
K-IS-NSK 
KA— I -MASK 

FO— 9M!K>— SUB  ' '* 

IF ! HMAX  .NE. OSCMAX  .AND. SLOPE. £0. 1 .£30) 

*■  P(KA)-PS*£XP!HH/2.-(SGRT(FI  >*SQRT!FO  >)*0.2S*ALOG!FI/F0 ) ) 
IFiHMAX.EO.OSCMAX.OR.SLOPE.NE. I .£30)  PHCA)-PS*SLOPE*iZM!IO-ZM< TS) > 
IF(P(KA)*P{  I )XT.O..OR.P!KA).EQ.O.)  NOOES-NOOES-I 

‘ *»‘r  > »<  TPT-ZM 1 K > ) .GE.O.)  HMAX— AM  IN  I (HMAX.  OSCMAX) 
PINT-HH*0.3*PS**2 
HTST-AMINI (HTST, HMAX) 

Gtt  TO  13a  - ..  ..  ;..  • . 

K-K-NSK 

ka-ka*nask  ..  . .... 

IFIKA.GT.NP)  Oa  TO  700  ' 

FO— MMiO—SUB  ....  „ ...  „ _ . . ■ 

IFiABS! I .— FF*FO) .LE.O(VCX)  GO  TO  400 

P(KA)— <P(KA-NASK3*<2.*IO.*FF*Fl  >— ?OCA-2*NASO*!  I w—FF*F2))/ 

*■  ! I .-FF*FO) 

kmax-ka 

IF(P(XA)*P(  KA-NA3K)  .LT.O..OR.  ?<XA ) .£0.0. ) NOOES-NOOES-I 
IF!  (ZM(K*N5IC)-rPT>*(  TPT-ZM(X)  J.GE.O.)  HMAX— AMIN  I (HMAX.!)SCMAX) 
IF!A8SIFa).LE.0IVCX)  HTST-HMAX 

>MXAMMXJQV.AX.nrynfl  MTtT.tHIIII  mi  — , 


140 


I 


it.. 


I'j 


- J 


. . 2 


fl 

h 


KT*. 


Page- 44 


* )),HMAX>  5 

(F<riH.GT.0.3*HTST)  CO  TO  ISO 

IF<KA.LS.2*NA5K.oa.lCA*2*NASIC.GT.NP.OR.iCA*NASIC.OT.IMID-H>  GO  TO  160 
lF<<ZM(IC*2*NSlC)-rPT)*<TPT-ZM(IC)).Lr.O.)  GO  TO  147 
HMAX*AMINICHMAX,0$CMAX) 

HT5f*AMINI (HTST.HMAX)  ~ • 

IF<HH.GT.0.5*HTST)  GO  TO  130 
147  PINT*PINT*I .5*«H*P<KA)**2 
HH*HH*2. 

FP«FF*4. 

NSK-NSX*2 

NASK*NASIC*2  - - - - 

C aRITEtd, 164JNASK 
C 164  FORMAT<IX.-NASlC-*,I5,*  67*)  - 

GO  TO  170 

150  IF<HH.LE.HTST.OB.ICAr-lM.I0.GE.4)  GOTO  160 
P INT»P I NT*.  5*tW*P  < ICA  > **2 
152  IF<NASK.LT.2)  GO  TO  300 
HH-HH/2. 

FF-FF/4. 

NSK-NSK/2 

NASK-NAS1C/2 

F2*VU(IC-NSIC)-SUB 

IFU8S12.*I0.*FF*F2).L£.0IVC!C)  GO  TO  600 
PUA-NAS1C)«(P<ICA)*<  l.-FF*F0>*P<KA-2*NASK)*U..-FF*FI  ))/ 

♦ U.*I0.*FF*F2> 

IF(HH.LE.HTST)  GO  TO  153  i 

FI-F2 

- IF1NASIC.LT. 2.  )NRITE(6, 1631NASIC  » { 

163  F0RMATIU.*NASK-*,13,*  81*) 

GO  TO  132 

155  PlNT»PINT*.5*t«*P<ZA)**2  . 

oo  to  1 7a 

160  F2-FI  ~ i-  ~ n . - 

PUtT»PIMT*Wi*POCA)**2  - ’ 

170  FI-FO - - —x*.  :;-.r  — 

IFOCA— IMID.LT.*)  GO  TO  140 

■ PIMT*PINT-0.S*HH*PtlCA)*»2  - .~  l 

return  5 

500  mRITE(6»  1300)  HTST.IC.ZMIIO  -*  - • - -•  - ~ 

1500  F0RMATl*MAA  STPSIZ-*.P14.6.*  TOO  SHALL.  R-*.I6»*  Z**.F8.2) 

STOP  • — • . . — r • : • ! 

600  i«RITS(6.I600)  K.ZMCIC) 

1600  FOBMATC«0IV  BY  0 At  Z-*.F8.2)  ... 

STOP 

700  4RITEC6.L70Q)  K - 

1700  FORMAT! *OVERFLD»ED  STORAGE.  IC**.I6) 

end  ‘ "r  " . "■ 

SUBROUTINE  INTERP1IL.IU.P.NPI«.IC0FP>  • - — . 

DIMENSION  PCI  ) . ..  . 

REAL  ICO.ICOSO  ’ - ‘ ' - 

COMNON/GEN/UG.ZM<800l)^a^:ACT.0eLZ,lCO.0rVCr.S08.SUB.If0Sa.VM<  80 

oo  bo  k»i.iu  ' ’ " "•  ’•  : W-.  '• . ~ 

lF<PtIL-K*l>.N6.I.S60>  GOTO  8*  ,\r. .*_  * T .r: 

80  CONTINUE  - i --  . .:«•  • - . - 

GO  TO  PE  • ; t r i- • .'  i.- - --wi-'rvrr.-rr  - A**:  • — r — ~ ■■ “.  v 

8*  1 1 * I LHC*  I v*  .......  # .■  . . 

• DO  85  JC*IU»NP)flL  •••  - - ' «>**£•**•.  ’■  -~-T. 

IF<P<K).N(.1.E60)  GO  TO  87  ‘ ' V . : J 

85  CONTINUE  Vv  ' ^s-v- ' ’ r 

00  TO  98  ••  . rrn.-.»-i».- 

57  I2*<  < ** . F J ^ f j ' 

ihkofp.eq.0)  kc»i  • • ••y-  - 'rV-7rT 


■3+ Jr**' 
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1 


r nsk*i  . 

v DO  90  KX-II.I2  f*  >vU 

IFINSK.GT. I .AfC.PdOC)  .ME. I .E50)  GO  TO  92 

ifipooo.so.i  .£<m>  nsk-nsk*i 

00  TO  90 

92  KL-KX-NSK  < .-o—  - - - • • . — - 

KlHKK-1 

FF*(DSLZ»NSIO  **2/ 12.  - — . - 

94  NSX*NSK/2 
FF-FF/4. 

00  93  K-KL.KU.NSK 

IFlP<K).NE.I.E60>  CO  TO  93  

, L*KC*K*KOFF 

<•  F2*2.*I0.*<VMCL>-SUB)*FF  - --  ' ** 

IF(A8S(F2).L£.0IVCX)  GO  TO  97  ' 

P(K)*<P<K*N5K)*U  .-FF*(VM(L*NSK)-SUB))  • - 

* ♦P(K-N3K2*U.-FF*<VM(L-NSIC>-5UB).))/F2 

93  CONTINUE  - - - 

. IFtNSK.GT. I ) 00  TO  94 

90  CONTINUE  ■ 

RETURN  " • • 

97  *RITE<6.I097)  K,ZM(K) 

- 1097  FORMAT(*0IV  CHECK  IN  INTERP  AT  K**,13.*  Z**.F8.2> 

STOP 

98  NRITSU. 1098)  K 

* 1098  FORMATUTOO  4 AN I BLANKS  FOR  INTERP,  K-*. 13) 

STOP 

ENO 

SUBROUTINE  OERIVCP.I.Ia.ISGN.DER) 

01  MENS ION  P( I ) 

C0MM0N/0EN/U0.ZM<800I  >,Z0,rACT.0ELZ.K0.0lVCJC.SQ8.SU8.K0Sa.VM<  80 
•01)  .HI2.HI6.0UM<7) 

I2*IS0N*ISGN 

AI-0.S*<PCI*I3GN)-P(l-ISGN)) 

BI«HI2*<<VM(IA*I  )-SU8)*P<I*ISGN)-<VMi IA-I >-SUB)*P( I-ISCN ) ) 

- - A2*0.5*tPU*l2>-PU-I2»  — 

B2*HI2*<<VM<lA*2>-SUB)*P<l*I2)-<VMCIA-2>-SUa)*P(  I-I2)) 

• DER*HIA*<— A I *t.JS625*A2— 7.4*8  i— .425*82) 

RETURN 

-•  ENO  ' > .iV  : 

SUBROUTINE  TURN(ZLL«ZRR,ZT) 

CONNOH/OEN/ UC.ZIK  800 1 ) . ZO,  FAtT,  OEU  , KO,  0 1 VCK.  S08  , SUB,  KOSO. 

* VM<800J ),HI2. HIB.EAPMAX.OSC MAX. BC.TPT.SLOPE.EPS, B 

REAL  ICO. ICOSO  . 

OATA  ZMEET/ 337.  3323427/ 

El  (Z>-)./<  1 .006106525*1 .6774)  6-5*Z) 

VI(Z)*l.-El<Z)»*2 

E2(Z)-t  ,/t  I .*3.487656-5*2)  - - .. 

V2<ZJ*I.-E2(£)«*2 

ZL-ZLL  . . . - . 

ZR*ZRR 

IF<ZL.LE.ZMEET)0O  TO  25 
VL-VI(ZL) 

OO  TO  22  ••  -1  V-i  _ . 

25  VL-92CZL)  - 

22  IFCZR.LE.ZMEETIGO  TO  27  ...  \_  . 

VR-MI(ZR) 

OO  TO  28  

27  VR*V2<ZR> 

28  OO  10  I-l ,4a  . 

ZT-.5*<ZL*ZR)  — - . 

- I F<ZI. LEANEST)  GOT  TO  29-  

VT-MKZT)  - ..  . 

oo  to  3a  ••  - . - ....  . . ..  • 

- » VT*V2(ZT7 

30— LEU  RSI, VT*U0)  f.T.I  BN88UBI  BETtlRlL 


f 
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